Abstract. We present 3-d numerical magneto- 
hydrodynamic simulations of a buoyant, twisted magnetic 
flux rope embedded in a stratified, solar-like model 
convection zone. The flux rope is given an initial twist 
such that it neither kinks nor fragments during its ascent. 
Moreover, its magnetic energy content with respect to 
convection is chosen so that the flux rope retains its basic 
geometry while being deflected from a purely vertical 
ascent by convective flows. The simulations show that 
magnetic flux is advected away from the core of the flux 
T-H [ rope as it interacts with the convection. The results thus 
■ support the idea that the amount of toroidal flux stored 
' at or near the bottom of the solar convection zone may 
currently be underestimated. 
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1. Introduction 

The concept of buoyant magnetic flux tubes is an essen- 
tial part of the framework of current theories of dynamo 
action in stars, particularly in the case of cool dwarf stars 
such as the Sun. Results from studies of buoyant magnetic 
fl-ux tubes carried out within th e esse ntially 1-d thin flux 
tube approximation (e.g. Spruit 1981 and Moreno-Insertis 
198^ ) are consistent with the observed latitudes of emer- 



gence and tilt angles of bipolar regions on the surface of 
the Sun (Fan et al. |1994| and Cahgari et al. |1995D . More 
general 2-d simulations of flux tube cross-sections have 
shown that cylindrical tubes are disrupted by a magnetic 
Rayleigh- Taylor instability (e.g. Schiissler 1979|, T singanos 



198C, Cattaneo et al. 1990 



1995, Moreno- 



^ Matthews et al. 

Insertis & Emonet 1996| ). This renders them unlikely to 
reach the surface unless the presence of fieldline twist in- 
troduces a sufScient amount of magnetic tension to sup- 



press this effect (Emonet & Moreno-Insertis 1996, 1998 
and Dorch & Nordlund 1998[ ). On the one hand, 3-d simu- 



lations of buoyant, twisted flux ropes have confirmed sev- 
eral of the results from 2-d simulations (Matsumoto et 

1999| ), and have further shown 



al. 1998 and Dorch et al. 



that the S-shaped structure of a twisted flux tube as it 
emerges through the upper computational boundary is 
qualitatively similar to the sigmoidal structures observed 
in EUV and soft X-ray b y the Yohkoh and SoHO sate l- 
htes (e.g. Canfield et al. |l999| and Sterhng et al. |2000|) . 
Moreover, tightly packed (5-spots may be interpreted as 
the emer gence of highly twisted, kinking flux ropes (e.g. 
Fan et al. 1999). On the other hand, it has been suggested 
that the value of the critical degree of twist needed to pre- 
vent the Rayleigh- Taylor instability may be unrealistically 
high in the 2-d case, and a smaller twist may be sufficient 
in the case of sinusoidal 3-d magnetic flux loops (Abbett 
et al. ^000[ ). 

In the solar convection zone, buoyant flux structures 
are constantly interacting with the surrounding convec- 
tive downdrafts and updrafts, and the question remains 
whether the quasi-steady behavior that the flux ropes 



reach in the later phase of their rise in 2-d simulations 
(Emo net & Moreno-Insertis 199^ and Dorch & Nordlund 
1998 ) is stable towards perturbations from the surround- 
ings, and whether the results found for 3-d flux ropes mov- 
ing in a essentially 1-d static stratiflcation are valid in the 
more realistic case. 

In this paper, we present our first results regarding 
the behavior of buoyant, twisted flux ropes embedded in 
a fully dynamic model of solar-like convection. 



2. Numerical model 

The set-up of the model is twofold, consisting of a snapshot 
of a time-dependent, but statistically relaxed "local box" 
convection zone model (sandwiched between two stable 
layers), and of an idealized twisted magnetic flux rope. We 
solved the full resistive and compressible MHD-equations 
on a staggered mesh of 150 vertical x 105^^ horizontal grid 
points, using the method by Galsgaard and others (e.g. 
Galsgaard & Nordlund [1997| , Nordlund et al. 
dp 
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u, P, B, and e represent the density, velocity, 
pressure, magnetic field, and internal energy respectively. 
rj and t denote the magnetic diffusivity and the viscous 
stress tensor; the source terms Qvisc, Qjouic, and Qrad refer 
to the viscous. Joule, and diffusive heating. In the upper 
part of the domain, Qrad includes an additional term that 
provides for a simple, isothermal cooling layer. 

The computational method employs a finite difference 
staggered mesh with 6th order derivative operators, 5th 
order centering operators, and a 3rd order time-stepping 
routine. The diffusive terms are quenched in regions with 
smooth variations, to reduce the diffusion of well-resolved 
structures. Magnetic Reynolds numbers in non-smooth re- 
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Fig. 1. The vertical velocity in slices of the model convection 
zone in the initial state at two different depths: at the surface 
(left) and at the depth of the flux rope (right). 



gions are of the order a few times 10^, but can be much 
higher in smooth regions. 

The computational box is horizontally periodic with 
sides of 250 Mm and a height of 313 Mm (of which 166 
Mm is the convection zone, that covers 2.5 orders of mag- 
nitude in pressure). The flows are turbulent through-out 
the convection zone, and the kinetic energy spectrum dis- 
plays a power law at intermediate wavenumbers (fc w 3 - 
10). As it is typical for over-turning stratified convection, 
a cellular granulation pattern is generated on the surface 
of the convection zone (Fig. |l|). The typical length scale 
of this pattern is about 50 Mm, somewhat larger than 
the canonical size of 32 Mm of solar super-granules (e.g. 
Leighton et al. 1962). The typical velocity of the granula- 
tion is 200 m/s in the narrow downdrafts at the surface 
and slightly less in the upwelling regions, close to what is 
found for solar super-granulation (e.g. Worden & Simon 



1976) 



We choose an initially isentropic flux rope with a buoy- 
ancy of 1/7 (with 7 = 5/3) lower than the case of tempera- 
ture balance, where the buoyancy is 1//3 (with (3 being the 
classical plasma beta). This is computationally advanta- 
geous, since we avoid the costly process of perturbing the 
flux rope from a state of mechanical equilibrium. 

The initial twist of the flux rope is given by 

= Boe-('-/^)' and = a(r/R)B„ (5) 
where B^ is the parallel and B^ the transverse component 
of the magnetic field with respect to the rope's horizon- 
tal main axis. The coordinate system is chosen so that z 
corresponds to the initial axial direction of the rope. 

The wavelength A of the flux rope is equal to the hor- 
izontal size of the domain so that A = 3.2 Hpo at the 
initial position of the rope. Thus, the flux rope is not un- 
dular Parker-unstable even though the stratification per- 
mits this instability for longer wavelengths (Spruit & van 
Ballegooijen |1982 ). The rope is initially twisted, with a 
pitch angle (at r = R) of ipYi — arctan(Q;) and a radius 
Ro = 0.177 Hpo which corresponds to a half- width at half- 
maximum of Bz (henceforth HWHM) of ~ 0.1 Hpq. 

To avoid problems associated with the large ratio of 
thermal to dynamic time scales, our convection model has 



a much higher luminosity than the Sun, and thus, all vari- 
ables must be scaled to compare with solar values. The 
choice of the magnetic field strength is somewhat prob- 
lematic in this regard. The ratio of kinetic to thermal en- 
ergy density ck/g is much larger in our model than in the 
Sun (though the convective flows remain subsonic with an 
average Mach number of 0.01). This requires a choice of 
[3 that is smaller than its presumed value at the base of 
the solar convection zone so that the ratio of magnetic to 
kinetic energy density cm / ck is the proper order of mag- 
nitude. However, a small (3 is what is needed so that the 
time it takes to complete a simulation is not prohibitively 
long. We choose j3 = 100, which yields a solar-like cm/gk 
of 100. 



3. Results 

We have performed several fully convective 3-d simula- 
tions, as well as a number of 2-d convection-less simula- 
tions. The results of the latter agree with previous 2-d find- 
ings, and are used for reference in the following. We dis- 
cuss results from a 3-d simulation with -0^ = 45° (a = 1). 
In that case, the degree of twist is small enough to pre- 
vent the onset of the kink instability (the linear growth 
rate vanishes for a = 1, e.g. Fan et al. 1999), yet it is 
large enough to prevent the onset of the Rayleigh- Taylor 
instability. Thus, the rope retains its cohesion without dis- 
torting its shape by any of these two instabilities, and we 
can focus our attention on the effects of the convective 
flows on the rope. 

Fig. ^ compares our 3-d simulation to a correspond- 
ing 2-d convection-less reference simulation, and a simple 
analytic flux tube. As the 3-d rope rises, convective flows 
perturb its motion, preventing it from entering a well- 
defined terminal rise phase with a constant rise speed, as 
in the 2-d reference simulation (see Fig. ||a) . The rope re- 
mains straight and the maximum excursion of its axis, at 
the end of the simulation, is ~ 0.04 A, where we define the 
rope's axis as the set of positions along the rope, where 
the magnetic field strength is maximum. With the chosen 
super-equipartition axial field strength, the main action of 
the large-scale convective flows is to push the rope both 
left and right of the central plane (Fig. see also the 
mpeg- movie at Dorch [2001 ), while the effect of the small- 
scale downdrafts (of the order of the rope's radius) is to 
locally deform its equipartition boundary. 

Initially the rope is located in a general updraft region. 
This explains why the rise speed of the rope is slightly 
greater than that of the 2-d reference simulation, which 
reaches a terminal speed of ~ 0.1 wao (i'ao being the 
Alfven speed at the initial position of the rope). Never- 
theless, the two ropes arrive at the same final height at 
the end of the runs (though in the 3-d case, we note that 
the rope follows a longer path) . Our 3-d rope also expands 
more quickly than the rope in the 2-d simulation, and its 
rate of expansion is closer to what is expected from an adi- 
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Fig. 2. a) height of the flux rope as a function of time (stars). Also plotted are the corresponding results from a 2-d reference 
simulation (diamonds). The straight line corresponds to the average speed (0.1 vao) in the rise phase, b) drift of the flux rope 
in the meridional plane, c) expansion of the flux rope (stars) with the corresponding result from the 2-d reference simulation 
over plotted (diamonds). Also plotted is an analytical expression (solid line, see text), d) magnetic field strength as a function 
of height. 



abatically expanding, non-stretching tube with constant 
flux {B/p =const.): 

-l/2Va 



R(x) = Ro 1 - Va 



H 



PO 



(6) 



Fig. 1^ shows the rope's characteristic size Rhwhm defined 
by the average between the vertical and horizontal HWHM 
along its axis (the short period oscillations due to differ- 
ential buoyancy, that are not well-resolved, have been fil- 
tered out by smoothing over two grid points). As the rope 
rises and expands, its magnetic field strength Be, here de- 
fined as the average axial field strength along the rope, 
decreases at a rate close to that determined by Eq. (^. 
At later times, the field strength of the 3-d and 2-d ropes 
decrease at nearly the same rate (Fig. |2|i). The deviation 
can be attributed to the fact that, during its ascent, a sig- 
nificant amount of the magnetic fiux within the 3-d rope is 
lost to its surroundings. This is illustrated by Fig. ^ (left), 
which shows the total normalized magnetic flux within 
the HWHM-boundary $i as a function of time for both 
the 2-d and 3-d ropes. We note that as the 3-d simulation 
progresses, the total flux-loss from the computational do- 
main is only 0.3%. The flux content of the rope, however, 
decreases much more quickly. 

Also shown in Fig. || (right) is the magnetic flux exter- 
nal $e to the rope both above and below its center and 
respectively. Since the sum ^e + ^i is nearly conserved, 
as $i decreases, = + must increase by an equal 
amount. However, the distribution of the flux-loss is not 
symmetric: more flux is lost to the surroundings below the 
rope than above. The e-folding time of the increase of flux 
in the lower domain is ^ 20 ta (with ta = Rq/wao)- 
This asymmetry also exists in the 2-d reference simula- 
tion, even though the total flux-loss is much smaller in that 
case. The asymmetry is a result of two factors. First, as 
the rope rises, the total volume above it decreases, while 
the volume below it increases. Second, there is an anti- 
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Fig. 3. Left: magnetic flux within the rope $i (stars), the 
corresponding quantity in the 2-d simulation (diamonds) and 
an analytic flt (solid curve), Right: the normalized flux outside 
and above the center of the 3-d rope $u (stars) , and below, $; 
(triangles) . The same quantities are shown for the 2-d reference 
simulation (solid and dashed curves respectively). 



symmetry of the relative velocity across the rope. When 
the rope ascends, there is a tendency for flux to be ad- 
vected towards the rope near its apex, and transported 
away from the rope in its wake. The more pronounced 
asymmetry in the 3-d case can be attributed to the pump- 
ing effect tha t tran sports the weak fi eld do wnwards (Dorch 
& Nordlund ^00l| and Tobias et al. ^00l[ ). 



We have defined the flux rope as the magnetic struc- 
ture that lies within the HWHM-boundary. This boundary 
is not, however, a contour that moves with the fluid in the 
classical sense: the flux within the latter kind of contour is 
naturally conserved (neglecting resistive effects) and equal 
to the total flux $o = ^'e + ^'i- The HWHM-boundary is a 
convenient way of defining the flux rope and a characteris- 
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tic size Khwhm, that behaves more or less as it is expected 
from the analytical expression Eq. (|^). The evolution of 
the flux within the rope's core $i is determined by 

$i = - i" Av X B • dl, (7) 

./boundary 

where Av is the difference between the fluid velocity and 
the motion of the HWHM-boundary. The average "slip" 
Av of the rope's boundary in the simulation is only a small 
fraction of the rise speed, and varies between the range of 
plus or minus a few times 10~^ and 10~^ z;ao- 

Making the rather crude assumptions that the bound- 
ary only moves radially relative to the fluid and that the 
circumference of the boundary is circular (which it is not), 
Eq. (0) reduces to 

= -TrRhwhm Av Be, (8) 

where Be is the field strength at the center of the flux 
rope. Integrating Eq. (||) numerically with Khwhm and Be 
determined from the simulation (see Fig. ^), the result 
is a perhaps surprisingly good fit to the actual flux- loss, 
see Fig. || (left), if Aw is set to 3 10~^ wao throughout 
the time span of the simulation except for a short interval 
of ~ 3 TA around t = 30 ta, where Av = — 1G~^ vaq, 
when the rope passes from one updraft to another (see 
the discussion below). 

4. Discussion and conclusions 

The 3-d numerical simulations show that the interaction 
of a buoyant twisted flux rope with stratified convection 
leads to a considerable loss of magnetic flux from the 
core of the flux rope (as defined by the rope's HWHM- 
boundary) . 

The initial position of a flux rope in the convection 
zone is significant for the subsequent detailed history of 
its rise: with the present convective flows and the initial 
location of the flux rope, most of the rope starts out lo- 
cated inside or close to a convective updraft. Thus, the 
ascent of the rope is likely to be influenced by this fact, 
and we are therefore not able to draw any conclusions on 
the detailed path of its rise. However, in the course of the 
simulation, the flux rope rises 96 Mm, and loses about 25% 
of its original flux content. This, ceteris paribus, leads to 
an increase in the amount of toroidal flux that must be 
stored at the bottom of the convection zone during the 
course of the solar cycle. 

In the Sun, toroidal flux ropes rise about 200 Mm 
through the convection zone before emerging as bipolar 
active regions. One may thus expect them to lose even 
more of their initial flux, which would then be pumped 
back down toward the bottom of the convection zone. We 
can quantify this subsequent flux-loss by assuming that 
Eq. (^ is valid through-out the rise, that the ropes ex- 
pand according to the simple analytical estimate of Eq. 
(^, and that the ratio of the slip Av to the rise speed 
remains constant. Given these assumptions, the flux-loss 
at a height of 200 Mm is 26 % of the initial flux, i.e. not 



much more than in our simulation. However, the relative 
slip may not remain constant throughout the rope's rise. 
For example, Av and thus $i changes at the time around 
i = 30 TA, which corresponds to the time when the rope 
is at its maximum (rightward) excursion from a vertical 
ascent (see Fig. ||b). At that time, the rope exits the con- 
vective updraft with which it was initially associated, and 
enters a different ascending "plume" to the left of the its 
original position. This leads to a transient compression of 
the rope {Av < 0, in the simplifled expression Eq. ||). Af- 
ter entering the new plume, the average slip returns to its 
previous positive value for the remainder of the rise. 

Petrovay & Moreno-Insertis ( 1997 ) suggested that tur- 
bulent erosion of magnetic flux tubes may take place 
within the solar convection zone due to the "gnawing" of 
turbulent convection. They propose a mechanism whereby 
a flux tube is eroded by a thin current sheet that forms 
spontaneously within a diffusion time. That we do not see 
a loss of flux via this type of enhanced diffusion should 
not, however, be taken as a dis-proof of the feasibility 
of turbulent erosion: it requires the turbulence to be re- 
solved down to much smaller scales £ <C A, than in our 
simulations. Instead, the flux-loss is completely due to the 
advection of flux away from the core of the flux rope by 
convective motions. Most of the flux that is "gnawed-off" 
ends up in the trailing wake and some of this flux is mixed 
back into the upper layers by ascending flows. We specu- 
late that both types of flux-loss may take place simultane- 
ously in the Sun, and as a result, the amount of toroidal 
flux stored near the bottom of the solar convection zone 
may currently be underestimated. 
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